Relevance of Abelian Symmetry and Stochasticity in Directed Sandpiles 
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We provide a comprehensive view on the role ol Abelian symmetry and stochasticity in the 
universality class of directed sandpile models, in context of the underlying spatial correlations of 
metastable patterns and scars. It is argued that the relevance of Abelian symmetry may depend on 
whether the dynamic rule is stochastic or deterministic, by means of the interaction of metastable 
patterns and avalanche flow. Based on the new scaling relations, we conjecture critical exponents 
for avalanche, which is confirmed reasonably well in large-scale numerical simulations. 
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Since a prototype of sandpile models was first intro- 
duced by Bak, Tang, and Wiesenfeld, lots of its variants 
have been tested and become successful over the last two 
decades in figuring out the underlying common mecha- 
nism of ubiquitous scale invariance in nature [l| . In such 
models, grains are slowly added, redistributed (toppled) 
instantly whenever the instability threshold is overcome, 
and finally dissipated at boundaries. It has been the most 
interesting issue and under debate whether the universal- 
ity class of critical avalanche dynamics can be changed 
by the modification of local toppling rules, such as the 
breaking of Abelian symmetry [2J and the consideration 
of stochasticity [3J|, with still conflicting numerical re- 
sults |4j . Abelian symmetry here means that the order of 
toppling the unstable sites does not affect the final state. 
Contrary to undirected models, directed sandpile models 
(DSMs) with a preferred direction of toppling turns out 
to be more tractable analytically as long as they have 
Abelian symmetry [fj|, |g, [7| . It is because metastable pat- 
terns in the Abelian DSMs are fully uncorrclatcd. 

Once the Abelian symmetry is broken in DSMs by 
some specific way, long-range spatial correlations emerge 
in their metastable patterns. Such correlations are often 
observed in nature, like a fractal structure in the crust 
of the earth formed by seismic events. Two non- Abelian 
DSMs with spatially correlated metastable patterns were 
introduced by Hughes and Paczuski [8[ for a stochastic 
version and by Pan et al. [9( for a deterministic version. 
In the stochastic version it is claimed that Abelian sym- 
metry is not relevant to avalanche dynamics, while in the 
deterministic version it is. Although this difference might 
be attributed to the existence of stochasticity, it is not 
clear enough to say which factor governs the scaling prop- 
erty of metastable patterns. Therefore, it is quite crucial 
to clarify the role of spatially correlated metastable pat- 
terns in the universality class of DSMs, which has been 
hardly discussed up to now. 

In this Letter, we discuss how the critical avalanche 
dynamics of non- Abelian models are entangled with spa- 
tially correlated metastable patterns. Based on the 
formation of metastable patterns and scars (trace of 



avalanche boundary sites) with the mapping onto particle 
dynamics, we give intuitive arguments about the scaling 
relations in terms of scar exponent, and conjecture a pos- 
sible scenario for the universality class in DSMs. Finally, 
we reinterpret the earlier known results for the Abelian 
case by our conjecture, and confirm those for non- Abelian 
case by large-scale numerical simulations with various 
data analysis techniques developed so far. 

Consider DSMs defined on a (1 + l)-dimensional tilted 
square lattice of size (L,T). The preferred direction 
of avalanche propagation is denoted by the 'layer' t = 
0, • • • , T — 1 with open boundary conditions, and the 
transverse direction by i = 0, • • • , L — 1 with periodic 
boundary conditions. Initially, to each site of the lattice 
an integer value (the number of grains), Zi(t) g [0, z c ), is 
assigned, where we set the instability threshold z c = 2. 
Given a stable configuration where all sites are stable, 
new grains are added one by one at a randomly cho- 
sen site on the top layer, Zi(0) — > z,-(0) + 1, until one 
of them becomes unstable. For any unstable site with 
Zi(t) > 2, grains at that site topple down to its left and 
right nearest-neighboring sites on the next layer, t + 1: 



Zi(t) -> Zi(t) - Aii, 
Zi±i{t + 1) — > z l± i(t + 1) + A M ±i, 



(1) 



where A;, = A; : ,_i + A^j+i (the local conservation of 
grains). Toppled grains at the unstable sites on the bot- 
tom layer t = T — 1 are dissipated out of the system. 
Only after another stable configuration is recovered by a 
series of toppling events, denoting an avalanche, a new 
grain is added to keep generating another avalanche. By 
setting {Ay} one may consider several variants of DSMs. 
In contrast to Abelian DSMs where An is constant (used 
to be set as z c ), we set An = z,-(£) as non- Abelian DSMs. 
All grains at the unstable site topple to the next layer 
and the toppled site becomes completely empty. For any 
given Ai», the values of A^±i can be determined in cither 
stochastic or deterministic way. Besides the well-known 
Abelian deterministic or stochastic DSMs [a, |6[ (AD/AS 
in short) and the non-Abelian stochastic DSM |8|] (NS 
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FIG. 1: Toppling rules of non-Abelian deterministic DSMs 
on a (1 + l)-dimensional tilted lattice: (a) aND, (b) bND, 
and (c) cND. The gray colored grains and arrow represent 
the state before toppling and the black colored ones represent 
the state after toppling, respectively. 



in short), we explore the following three versions as the 
non-Abelian deterministic DSM (ND in short): 



(i) A M± i 
(ii) A, 



k if Zi(t) = 2k 

fc + £t±i,<n(t) H Zi(t) = 2k + 1, 

k iiz t (t) = 2k 

i k + 5 l±1 , i+1 if Zi (t) = 2k + 1, 

(iii) Ai,i±l = : Zi{t)/2, 



where A; is a positive integer and Sij denotes the Kro- 
ncckcr delta-function. Wc call (i) the alternatively biased 
version (aND) [9j, (ii) the fully biased version (bND), and 
(iii) the continuous version without bias (cND), respec- 
tively. For the aND, an 'arrow' Oj(i) of each site initially 
points to one of its neighbors, say i — 1 in Fig. Q] (a). 
Whenever each grain is toppled at that site, the direc- 
tion of the arrow flips to the other neighbor, see Fig. [TJ 
which shows the case of Zi(i) = 3. 

Each avalanche can be characterized by the following 
quantities: mass s (the number of toppled grains), dura- 
tion t (the number of affected layers), area a (the num- 
ber of distinct toppled sites), width w (the mean dis- 
tance between left and right boundaries of avalanche), 
and height h (the mean number of toppled grains per 
toppled site). The avalanche distribution functions in 
DSMs show no characteristic scale except for T as long 
as L is sufficiently larger than the maximum width. They 
follow the simple scaling form as P(x) ~ x~ Tx f (i/T Dl ) 
for x € {s,t,a,w,h}. Moreover, two quantities x and 
y scale as (y) ~ x ly:c with j yx = 



D„ 



from 

P(x)dx — P(y)dy. Taking full advantage of the rela- 
tions, D t = 1 and (s) ~ T in DSMs, with the reasonable 
assumption of compactness of avalanche, i.e. a ~ wt 
and s ~ ah, we obtain the following scaling relations: 
Ixt = D x for any x, D s (2 - r s ) = 1, D a = D w + 1, 
and D s = D a + Dh- As a result, there are only two 
independent exponents left. Concerning the metastable 
state, we define two scaling exponents more. Along the 
propagation direction, one can measure grain density as 
pit) = (j; J2i z i(t)) ~ t~ a w ith the grain density expo- 
nent a for large t. The scar density yO sc (i) and the scar 
exponent a sc arc defined by the same definition of grain 



density only with Zi(t) replaced by hit), which takes a 
value of 1 for trace or otherwise. The scar exponent 
is immediately related to the avalanche width exponent 
as a sc = D w because the density of avalanche boundary 
sites is inversely proportional to the typical avalanche 
width, i.e. p sc (t) — w{t) . Since grains can remain 
only at the avalanche boundary sites for the non-Abelian 
DSMs, it is found that a = a sc . 

We give intuitive arguments on the interplay between 
avalanche flow and metastable patterns/scars in DSMs. 
Let us define N(t) as the number of grains transferred 
from the layer t to the next layer t+l within an avalanche, 
scaling as N(t) ~ w(t)h(t) ~ t D ~ +Dh . The evolution of 
N(t), avalanche flow, can be written as 



dN(t) 
dt 



N(t)-N(t-1)= J2 »i(*)> ( 2 ) 

iew(t) 



where rii(i) denotes the amount of the avalanche flow 
at each site (i,t), and the summation is over the sites 
between avalanche boundaries belonging to w. 

We begin with the Abelian case, for which the 
metastable patterns are fully uncorrelated. In the AD, 
it is well-known that Dh — by definition and D w = 
1/2 by mapping avalanche boundaries onto the random 
walks [5[ . The avalanche flow of the AD can be written as 
^- w r\. An uncorrelated noise r\ of zero mean and unit 
variance denotes the fluctuation of grain density. In the 
AS with the same 77, the bulk contribution to avalanche 
flow plays a crucial role, so that we get ^£- w \pwr\ [0, Lflj- 



Here y/w represents the fluctuation of the number of top- 
pled sites in w when the topplings are uncorrelated. The 
lack of correlation in metastable patterns again leads to 
D w = 1/2, so Dh = 1 ~® W = I/ 4 - In relation to the 
scar density, we get a sc = D w = 1/2 but a — for the 
Abelian case. It is found that the measurement of a sc is 
more efficient than that of D w . From now on, we suggest 
to measure a sc as one independent exponent in DSMs. 

For the non-Abelian case with the spatial correlation of 
metastable patterns or scars, the fluctuation of grain den- 
sity is numerically found to scale as pit). Furthermore, 
p sc plays the same role as p, i.e. a sc = a. In the NS [8|, 
keeping y/w due to the stochastic nature of toppling just 
as in the AS, we can write —nr ~ \fwp and get the new 



scaling relation, Dh = 1 — 



In contrast to the NS, 



in the ND, all sites of w contribute to avalanche flow so 
that ^n- f=a wp. We get the different relation Dh = 1 — a. 
Therefore, all avalanche exponents in the non-Abelian 
DSMs can be obtained from the grain density exponent 
a (or a sc ) as shown in Tabic [T] 

The new scaling relations with a enable us to clarify 
the effect of metastable pattern with a w 0.45 in the NS 
on its avalanche dynamics. Interestingly, if a = 1/2 is 
assumed with systemic errors and/or possible logarith- 
mic corrections, the NS has exactly the same avalanche 
exponents as those of AS: t s = b~) = 10/7 an d n = 




FIG. 2: (color online) Typical metastable patterns for non- 
Abelian case: The occupied sites are shown as black dots and 
the typical shapes of dissipative avalanches consisting of the 
toppled sites as blue (or gray)-shaded areas on a lattice with 
L = 150 and T = 250. Here (a) aND, (b) bND, (c) cND, and 
(d) NS, respectively. 



2 — ^ — 7/4. Another scenario for a = 1/2 can be found 
by mapping metastable patterns onto the space-time con- 
figuration of 2 A — > A coagulation-diffusion model defined 
in d = 1, where the particle density decays as t -1 / 2 [101 ] . 
One can say that the NS belongs to the same universality 
class as the AS in the following sense: For Abclian cases 
the flowing avalanche can sweep and lose many grains at 
the same time due to the uniform grain density. On the 
other hand, for non-Abelian case the flowing avalanche 
can sweep only a few grains due to the power-law decay- 
ing grain density and leave few grains behind by taking 
all grains at the toppling sites. In other words, the scal- 
ing property of N(t) is apparently unaffected by the grain 
density as long as the toppling rule is stochastic. 

From Dh = 1 — a in the ND, we find that the avalanche 
exponents for mass and duration have the mean-field 
(MF) values, independent of a, i.e. r s = 3/2 and r t = 2, 
whereas other exponents depend on a. We point out 
that one should consider other avalanche exponents as 
well as those of mass and duration in order to discuss 
the universality class of the ND. If a = D w = 1 is as- 
sumed from the linear behavior of avalanche boundaries 
(scars) as shown in Fig. [5J all avalanche exponents turn 



to the MF values, except for the case of width ll| . This 
may also correspond to the MF behavior of coagulation- 
diffusion model, where the particle density decays as t~ l 
in d > d u = 2 [10J. This peculiar 'MF' behavior of three 
ND versions appeared in the low dimensional system can 
be understood by considering the shape of N(t) with 
its width and height. We now focus on how avalanche 
boundaries behave linearly, which implies D w = 1. The 
ND toppling rules we considered suppress the fluctua- 
tions of height profile of the flowing avalanche more than 
the stochastic one does, which leads to spread grains 
wider and make the avalanche boundaries grow faster, al- 
most ballistically. This positive feedback enables D w = 1 
to be larger than 1/2 for all other DSMs. Moreover, we 
like to note that the resultant Dh = indicates the MF 
behavior for the non-Abelian case, whereas Dh = for 
any dimension in the AD. 



We performed extensive numerical simulations for all 
DSMs to confirm our conjecture about the avalanche ex- 
ponents in terms of the scar exponent, a sc = a, up to 
T = 2 13 and L = T/2 (T = 2 15 or L = T in some 
cases). We measure the avalanche exponent set {t Xi D x } 
of all x, a, and a sc using the moment analysis and the 
conventional successive slope techniques of avalanche dis- 
tributions, for about 10 9 avalanches at the steady state 
after the transient period. The spatially correlated scars 
are observed in all DSMs with the nonzero values of a sc , 
while the spatial correlations of metastable patterns are 
only observed in non-Abelian DSMs with the nonzero val- 
ues of a. To validate the stability of the scar exponent, 
we plotted the effective scar exponent, [a sc ]off , as a func- 
tion of 1/t for all DSMs. As shown in Fig. 3(a), there 
seems to be two asymptotic values, 1/2 and 1 in the large 
t limit, with quite long/unusaul initial transient behav- 
iors and finite-size corrections. For the non-Abelian case, 
we also checked the possibility of logarithmic corrections 
to scaling in p sc and p, both of which behave qualita- 
tively the same. Thus, we only show p in Fig. 3(b) as 
pit) ~ t~ a (lnt)^, where the existence of linear parts in 
curves represent logarithmic corrections. 

We finally discuss the relevance of Abelian symmetry 
in DSMs. Based on our results, Abelian symmetry turns 
out to be irrelevant to the stochastic version only when 
a = 1/2 in the NS. The breaking of Abelian symmetry 
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FIG. 3: (color online) (a) The effective scar exponent [a sc ] c ff 
as a function of 1/t, and (b) double-logarithmic plots of p(t)t a 
versus lni for logarithmic correction checks, where a — 1 
unless noted. For the bND, we cannot exclude the possibility 
of a = 2/3 either. Here the data of each model were obtained 
from 10 8 or more avalanches on a lattice of L = 2 14 and 
T = L, except for the cND where T = 2 13 . 



TABLE I: Avalanche exponents {t x ,D x }, grain density exponent a and scar exponent a sc in (1 + l)-dimensional DSMs with 
our conjecture. Note that a ac = a for the non-Abelian case. We check the value of K x = D x (t x — 1), which is a universal 
constant for any x [12|], and show its averaged value over x excluding K w and Kh due to their poor statistics. 

Model t s , D s 

Mean field (d u = 2) 3/2, 2 

Abelian 

Deterministic 4/3, 3/2 

Stochastic 10/7, 7/4 



Non-Abelian 
Stochastic 

- numerics 
Deterministic 

- aND 

- bND 

- cND 



T t , D t 


T a , D a 


Tw , Uw 


Th, D h 




(K*) 


ct, a sc 


2,1 


3/2,2 


3, 1/2 


oo, 




1 




3/2,1 
7/4,1 


4/3, 3/2 
3/2, 3/2 


2, 1/2 
5/2, 1/2 


oo, 
4,1/4 




1/2 
3/4 


0,1/2 
0, 1/2 


— - 1 

2 ' ± 


2fl+c0 ' x T u 


2±^, a 

2a ' 


4(1-0,) -, 

2-3c« ' x 


3a 

2 


1-f 


a 



2(3-a) 



1.43(1), 1.77(2) 1.78(1), 1.00(1) 1.53(1), 1.46(2) 2.74(2), 0.44(1) 3.18(2), 0.31(3) 0.77(1) 0.45(3) 



3/2,2 



2, 1 



2+a 



1 + Q 



- J —, a z — , 1 — a 



1.49(1), 1.97(3) 1.94(1), 1.00(1) 1.52(1), 1.88(2) 2.10(1), 0.87(1) 6.64(3), 0.07(1) 0.96(1) 0.86(3) 
1.43(1), 1.82(4) 1.79(1), 1.00(1) 1.48(1), 1.76(7) 2.10(1), 0.76(6) 5.90(9), 0.06(1) 0.81(2) 0.69(5) 
1.52(3), 1.99(4) 2.04(3), 1.00(1) 1.51(1), 1.95(3) 2.03(1), 0.86(2) 9.26(7), 0.06(4) 0.99(2) 0.91(11) 



in the deterministic version yields the MF behavior of 
avalanche dynamics even in a (l+l)-dimensional setup. 
In all NDs, we also confirm that the values of Dh are 
quite close to 0, which can be the sign of the MF be- 
havior as we argued. Furthermore, it turns out the NDs 
do not show any criticality in a (0+l)-dimcnsional setup. 
All numerical results arc listed in Tablc[I]with our conjec- 
ture. Only the results of the bND seem to be inconsistent 
with the values we conjectured. Such discrepancy may 
be attributed to the relevant effect of bias in toppling 
rule or relatively large logarithmic corrections. The va- 
lidity of our conjecture for other possible values of a or 
a sc is under investigation, with the role of toppling bias 
in DSMs " 
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In summary, we have explored the role of Abelian sym- 
metry and stochasticity in directed sandpiles, and con- 
jectured the new scaling relations for critical avalanche 
dynamics entangled with the underlying structure of 
metastable patterns and scars. Our conjecture provides 
clear guidelines on discussing the universality class in di- 
rected sandpile models. Moreover, our results provide es- 
sential information on analyzing the self-organized criti- 
cality in real systems as well as answering how ubiquitous 
long-range spatial correlations in nature can be developed 
and affect real avalanche dynamics. 
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